Self-similar approximations for trapped Bose-Einstein condensate 
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An approximate solution to the Gross-Pitaevskii equation for Bose-Einstein condensate in a spher- 
ical harmonic trap is suggested, which is valid in the whole interval of the coupling parameter, 
correctly interpolating between weak-coupling and strong-coupling limits. This solution is shown 
to be more accurate than the optimized Gaussian approximation as well as the Thomas-Fermi ap- 
proximation. The derivation of the solution is based on the self-similar approximation theory. The 
possibility of obtaining interpolation formulas in the case of nonspherical traps is discussed. 
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Bosc-Einstcin condensates of trapped atomic gases at low temperatures are described by the Gross-Pitaevskii 
equation (see reviews [1-3]). This nonlinear equation has no exact solution, and for practical applications one needs 
to resort to some approximations. The most commonly used such approximations are the Gaussian approximation 
and the Thomas-Fermi approximation. The former is asymptotically correct in the weak-coupling limit, while the 
latter, in the strong-coupling limit. Both these approximations are not accurate in the intermediate region of the 
coupling parameter. The Gaussian approximation can be improved by invoking an optimization procedure. Then its 
region of validity with respect to the coupling parameter essentially increases. However, it does not become exact 
in the limit of an infinite coupling. The Thomas-Fermi approximation, in addition to being incorrect at small and 
moderate coupling parameters, is always inadequate near the edge of an atomic cloud. There exist several suggestions 
[4-7] for correcting this approximation. But then the main advantage of the Thomas-Fermi approximation of being 
simple becomes lost. It would be desirable to have an approximation that would be valid for the whole range of the 
coupling parameters and would not require additional corrections. 

The aim of the present paper is to derive an approximate solution to the Gross-Pitaevskii equation with spherical 
symmetry, such that it would be accurate for arbitrary values of the coupling parameter and also would yield asymp- 
totically exact solutions in both limits of weak as well as strong coupling. We compare the derived solution with those 
obtained in an optimized Gaussian approximation and the Thomas-Fermi approximation and show that it is more 
accurate than the latter two. We also suggest the ways of generalizing the method for nonspherical traps. 

Let us consider a Bose-Einstein condensate of atoms with mass too in a harmonic spherical trap of frequency luq. 

It is convenient to work with dimcnsionlcss quantities, measuring the radial variable r = \Jr 2 + Ty + r 2 /lo in units 

of the oscillator length Iq = ^h/mQUa. The related dimensionless ground-state wave function ipo(f) depends on the 
variable r. The dimensionless coupling parameter g = ATra s N/lo is expressed through the ratio of the s-wave scattering 
length a s to the length l n , with N being the number of atoms in the trap. The radial wave function is defined as 
X(r) = V^r^o (r) . In this notation, the stationary Gross-Pitaevskii equation, for a harmonic trap, reduces to the 
eigenvalue problem 

HrX = E X , H r = -l^ + ^ + ^ x >. (1) 

In order to be self-consistent, the solution to the problem (1) has to satisfy the self-consistency conditions E = (%, H r \) 
and (x, x) = 1- To find an approximate solution to the problem (1), we shall employ the self- similar approximation 
theory [8-11], using the variant [12,13] designed for constructing crossover approximants satisfying the prescribed 
asymptotic conditions. Thus, we may notice that in the Hamiltonian there are two terms of different physical nature. 
The term r 2 /2 corresponds to the harmonic trapping potential, while the last term is due to atomic interactions. 
The contribution of these two terms is different for different values of the variable r. There exists a crossover radius 
r c = (g/271-) 1 / 4 , separating the axis r > onto two regions, where one of these terms plays the major role. For 
r <C r c , the interaction term prevails over the harmonic one, but for r 3> r c , the main contribution comes from 
the harmonic term. It is straightforward to check that for small r the solution to Eq. (1) reads as an expansion 
x(r) ~ ar + c 3 r 3 + c 5 r 5 , as r <C r c , while for large r it has the form \{r) ~ Cr exp(— r 2 /2), asr> r c . As is clear, 
the limit r — > is related to the strong-coupling limit g — > oo, when the interaction term dominates, while r — > oo 
corresponds to the weak-coupling limit g — > 0, where the harmonic term becomes prevailing. Our aim is to find a 
solution that would be accurate for all r and, respectively g, and would satisfy the prescribed asymptotic conditions. 

The idea of the method, we shall use, is to construct a self-similar interpolation between the asymptotically exact 
expressions. The name of the method comes from employing renormalization group in the form of group self-similarity 
[8-11]. All mathematical foundations of the theory and technical prescriptions are expounded in detail in Refs. [8-13]. 
The result of the self-similar interpolation between the asymptotic expressions is the self-similar approximant 

X* (r) — Cr exp ^— ^-^ cxp {ar 2 cxp(— or 2 )} . (2) 

Here C is a normalization constant and the parameters a and b arc defined by expanding Eq. (2) in powers of 
r and substituting this expansion into Eq. (1). Equating the coefficients at the powers of the same order gives 
a= 1/2 + (gC 2 -4nE)/12w, b= [2(1 - 2a) E- 2(1 - 2a) 2 - l]/20a. The energy E = E* in this approximation and the 
normalization constant C are defined by the self-consistency conditions, that is from the equations E* — (x*, H r x*) 
and (x*,x*) = 1. By this construction, the self-similar approximant (2) possesses the correct behaviour at small 
r <C r c and at large r 3> r c , thus, interpolating between asymptotically exact expressions. 

The accuracy of an approximation can be verified by calculating the local residual R(r) = (H r — E)x(r) and the 
integral deviation, a = [J °° {R^r^dr] 1 ^ 2 , defining, respectively, the deviation of an approximate solution from the 
exact one at each given point r and giving the integral measure of accuracy. Before calculating these, we shall recall 
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some other known approximations in order to compare their accuracy with that of the self-similar approximation 
and also to explicitly demonstrate what difference in such a physical quantity as the atomic density n(r) = X 2 ( r )/ r2 
results from the usage of different approximations. 

When the coupling parameter g is small, one may solve the nonlinear eigenproblem (1) by means of perturbation 
theory starting with the linear Hamiltonian. There exists an analytical continuation of linear modes to nonlinear 
stationary states [14]. Direct application of perturbation theory is valid only for asymptotically small coupling 
parameters g — > 0. In order to make perturbation theory relevant for finite values of g, it is necessary to invoke 
an optimization procedure. Thus, one comes to optimized perturbation theory. This approach was formulated [15] 
for treating the systems whose particles strongly interact with each other. The theory has been applied to various 
problems of quantum mechanics, statistical physics, condensed matter physics, and quantum field theory [3,11,15-19]. 
The optimization is realized by means of control functions [15]. In the present case, we may start from an approximate 
Hamiltonian with a harmonic potential u 2 r 2 /2 containing a trial parameter u. The ground-state eigenfunction of this 
Hamiltonian is of the Gaussian type, being 

/ 3\ V 4 

XG{r) = 2(^-\ rexp(-^r 2 ) . (3) 

Then, applying some variant of perturbation theory, e.g. Rayleigh-Schrodingcr theory, one can find a sequence 
{Ek(g,u)} of approximations for the energy. For instance, the first approximation is 

£l(9 ,„) = f(„ + I) + §„3/ 2 , ss _|_. w 

Control functions Uk(g) are defined so that to render the sequence {Ej c (g,Uk{g))} convergent. For example, the 
optimization condition dEk(g,u)/du = may be employed, resulting in the solution u = Uk(g), which for the case 
(4) gives the equation su 5 / 2 + u 2 — 1 = 0. If we stop at the first step of the optimized perturbation theory, then 
we get the optimized Gaussian approximation E G (g) = Ei(g,Ui(g)). Here we limit ourselves by this approximation, 
though higher-order corrections can also be obtained. Another popular approximate solution is the Thomas-Fermi 
approximation which is often used because of its simplicity. In that case, one ne glects the kinetic term in the Gross- 
Pitaevskii equation (1), which yields Xtf{t) = rQ{r 2 — r 2 )yj (2n/ g)(r 2 — r 2 ), where 0(-) is a unit-step function 
and ro = 2Etf- The energy is obtained from the normalization condition for the function, which gives Etf — 
i(15g/47r) 2 / 5 . Since the first of the self-consistency conditions is not satisfied, the Thomas-Fermi approximation is 
not self-consistent. 

To compare the accuracy of the approximations described above, we calculate the local residual, integral deviation, 
and the related energies by using, respectively, the self-similar approximant (2), optimized Gaussian approximation (3), 
and the Thomas-Fermi approximation xtf- The calculations show that the self-similar approximant (2) possesses the 
lowest residual, providing the most accurate solution for the Gross-Pitaevskii equation. In the case of the coupling 
parameter g = 25, the residual for the self-similar approximant gives < 0.01; for the optimized Gaussian 

approximation, < 0.4; while for the Thomas-Fermi approximation, it diverges near the edge of the atomic 

cloud. In the case g — 250, for the self-similar approximant, we have \R(r)\ < 0.4; for the optimized Gaussian 
approximation, |i?(r)| < 1.0; and the residual for the Thomas- Fermi approximation diverges near the classical turning 
point. The calculation of the integral deviation confirms that the self-similar approximant has the best accuracy. The 
deviation a for this approximant does not exceed 2.02 for all g > 0. It is very small, a <C 1, for the coupling g from 
weak to moderate. It slightly increases with increasing g, reaching a = 2.02 at g = 2411. Then it again decreases 
as g continues to increase. For instance, a = 1.04 at g = 2500. And again, a <C 1 for g — > oo. For the optimized 
Gaussian approximation, the deviation a slowly increases with rising g. Thus, a = 3.63 at g — 2500. Hence, this 
approximation works yet rather well even for sufficiently strong coupling parameters of the order g ~ 10 3 . In the case 
of the Thomas-Fermi approximation, the deviation a is infinite. 

Figure 1 presents the behaviour of the energy for the related approximants, with varying the coupling g. The 
optimized Gaussian approximation is asymptotically exact in the limit g — > 0, while the Thomas-Fermi approximation 
becomes asymptotically exact for g — > oo. The self-similar approximation interpolates between these two limits, being 
the best approximation uniformly valid for all g. The optimized Gaussian approximation is quite reasonable up to 
g ~ 10 3 . The Thomas- Fermi approximation is bad for low g and becomes reasonable after g ~ 100. Figure 2 illustrates 
the atomic density in different approximations as a function of the variable r for different couplings g. For g = 5, 
the self-similar and optimized Gaussian approximations are close to each other, being quite accurate solutions, which 
is confirmed by low residuals and small deviations a. The Thomas-Fermi approximation for such a relatively weak 
coupling is yet very inaccurate. For g = 25, the self-similar and optimized Gaussian approximations are still close to 
each other, and the Thomas-Fermi approximation starts approaching them. When g = 250, all three approximations 
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are of comparable accuracy, except that the Gaussian approximation slightly worsens near the trap center and the 
Thomas- Fermi approximation, at the edge of the atomic cloud. With varying the coupling from weak g — > to strong 
g — > oo, the self-similar approximation smoothly interpolates between the optimized Gaussian and Thomas- Fermi 
approximations, remaining the best approximation for all g, which is practically indistinguishable from the exact 
numerical solution of Eq. (1). A special caution is in order when the studied functions are small, as it happens 
for the asymptotic tails of our solution. Therefore, it is necessary to pay a particular attention to such asymptotic 
tails, checking if the relative difference is asymptotically small. Fortunately, for our case, the relative deviation of 
the self-similar approximation x*( r ) from the exact numerical solution tends to zero at large r as X*( r )/x( r ) — 1 — 
ar 2 exp(— br 2 ) — > 0, that is, the approximate and exact solutions asymptotically coincide. 

It is also worth commenting on the possibility of using the self-similar approximation technique for nonspherical 
traps. There are two ways of applying this method to equations containing more than one variable. One way is by 
reducing the equation in several variables to an effective equation of one variable, e.g. by means of an averaging 
procedure [3,14]. Another possibility could be by looking for an approximate solution to the three-dimensional Gross- 
Pitaevskii equation by constructing a trial function that is factorized with respect to its variables, as is done for an 
anisotropic trap in Rcf . [20] . It is also possible to invoke trial functions at the intermediate step for deriving analytical 
expressions for energies, which would interpolate between weak-coupling and strong-coupling limits. For this purpose, 
one should derive asymptotic expansions for the energy in the regions of small and large coupling parameters, and 
then construct self-similar approximations interpolating between these asymptotic expansions. As an example, we 
present the results of such an interpolation procedure for the ground-state energy of atoms in a cylindrical trap with 
the aspect ratio v = w z /w r , where ui z is the longitudinal trap frequency and u r is its transverse frequency. We give 
here the first three self-similar intcrpolative approximants for the energy measured in units of u> r , 



E{ = a (1 + AG) 2/b , 



E* 



a \(1 + A 1 G) 6/5 + A 2 G 2 



1/5 



£3 — «o 



{l + B l Gf /b + B 2 G 2 



11/10 



+ B 3 G S 



2/15 



(5) 



where G = 2gv/(2n) 3 / 2 , g = Ai:a s N/l r (l r = yJh/m a LO T ), and 

Aal /2 = 1.746928 , A ia f /6 = 2.533913(2 + v 2 f^ 



A 2 a% = 3.051758 



B ia 125/22 (2 + v 2 )° /m = 1.405455 (8 + 12^ 2 + //) 



.5/66 



5/6 



B 2 a 5/11 = 6.6 1 96 20 (2 + v 2 ) 



10/11 



B 3 al 5/2 = 5.331202 



ao = \ + v/2. 



Although the derivation of such formulas requires some work, but after being derived, they can strongly facilitate 
the overall consideration, since the energy is now given in sufficiently simply analytical form that is convenient to 
study with respet to the dependence on the coupling parameter g and the trap aspect ratio v. We have estimated 
the accuracy of these formulas in the intervals < g < 10000 and < v < 100. The maximal percentage errors are 
between 4 - 12% for E{- between 2 - 5% for E*: and of order 1% for E*. 

In conclusion, by employing the self-similar approximation theory [8-13], we have found an approximate solution 
of the Gross-Pitaevskii equation for a spherical trap. This solution, presented by the self-similar approximant (2), is 
compared with the optimized Gaussian approximation and Thomas-Fermi approximation. It is shown that among 
these three approximations the self-similar approximant provides the best accuracy for all couplings g, tending in the 
limits of weak and strong couplings to the corresponding asymptotically exact solutions. 



4 



[1] A.S. Parkins and D.F. Walls, Phys. Rep. 303, 1 (1998). 
[2] F. Dalfovo, et al., Rev. Mod. Phys. 71, 463 (1999). 

[3] P.W. Courteille, V.S. Bagnato, and V.I. Yukalov, Laser Phys. 11, 659 (2001). 
[4] F. Dalfovo, L. Pitaevskii, and S. Stringari, Phys. Rev. A 54, 4213 (1996). 
[5] E. Lundh, C.J. Pethick, and H. Smith, Phys. Rev. A 55, 2126 (1997). 
[6] A. L. Fetter and D.L. Feder, Phys. Rev. A 58, 3185 (1998). 
[7] P. Schuck and X. Vinas, Phys. Rev. A 61, 043603 (2000). 
[8] V.I. Yukalov, Phys. Rev. A 42, 3324 (1990). 
[9] V.I. Yukalov, J. Math. Phys. 32, 1235 (1991). 
[10] V.I. Yukalov, J. Math. Phys. 33, 3994 (1992). 

[11] V.I. Yukalov and E.P. Yukalova, Ann. Phys. (N.Y.) 277, 219 (1999). 

[12] V.I. Yukalov, E.P. Yukalova, and S. Gluzman, Phys. Rev. A 58, 96 (1998). 

[13] S. Gluzman and V.I. Yukalov, Phys. Rev. E 58, 4197 (1998). 

[14] Y.S. Kivshar, T.J. Alexander, and S.K. Turitsyn, Phys. Lett. A 278, 225 (2001). 

[15] V.I. Yukalov, Mosc. Univ. Phys. Bull. 31, 10 (1976). 

[16] W.E. Caswell, Ann. Phys. (N.Y.) 123, 153 (1979). 

[17] I. Halliday and P. Suranyi, Phys. Rev. D 21, 1529 (1980). 

[18] A.N. Sissakian and I.L. Solovtsov, Phys. Part. Nucl. 30, 1057 (1999). 

[19] M.B. Pinto and R.O. Ramos, Phys. Rev. D 61, 125016 (2000). 

[20] V.I. Yukalov, E.P. Yukalova, and V.S. Bagnato, Phys. Rev. A 56, 4845 (1997). 



•5 



Figure Captions 



Fig. 1. Condensate energy (in dimensionless units) as a function of the coupling parameter for the self-similar 
approximation (solid line), optimized Gaussian approximation (dotted line), and for the Thomas- Fermi approximations 
(dashed line). 



Fig. 2. Atomic density (in dimensionless units) as a function of the radial (dimensionless) variable, corresponding 
to the self-similar approximant (solid line), optimized Gaussian approximation (dotted line), and to the Thomas- Fermi 
approximation (dashed line) for different coupling parameters: (a) g = 5; (b) g = 25; (c) g = 250. 
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